/*
 * ElectrolyteSolution.cpp
 *
 *  Created on: 19 Aug 2011
 *      Author: Allan
 */

#include "ElectrolyteSolution.h"

// GeoReact includes
#include "Activity/AqueousActivity.h"
#include "Database.h"
#include "SpeciesMap.h"
#include "Utilities/Algorithms.h"
#include "Utilities/SpeciesUtils.h"

ElectrolyteSolution::ElectrolyteSolution()
{}

ElectrolyteSolution::ElectrolyteSolution(const vector<string>& species) :
species(species)
{
	// Determine the aqueous ions from the given in (species)
	ions = FilterIons(species);
	
	// Determine the aqueous complexes from the given in (species)
	complexes = FilterComplexes(species);
	
	// The species map object to determine the species indexes
	const SpeciesMap map(species);
	
	// Set the indexes of the ions
	idxions = map[ions];
	
	// Set the indexes of the complexes
	idxcomplexes = map[complexes];
	
	// Set the index of the water species
	iH2O = map["H2O"];
	
	// Determine the electrical charges of the ions
	zi = ElectricalCharges(ions);
	
	// Build the stoichiometric matrix between the aqueous complexes and ions
	vci = StoichiometricMatrixComplexesIons(complexes, ions);
}

const vector<string>& ElectrolyteSolution::GetSpecies() const
{
	return species;
}

const vector<string>& ElectrolyteSolution::GetIons() const
{
	return ions;
}

const vector<string>& ElectrolyteSolution::GetComplexes() const
{
	return complexes;
}

const VectorXd& ElectrolyteSolution::GetIonicCharges() const
{
	return zi;
}

const Index	ElectrolyteSolution::GetWaterIndex() const
{
	return iH2O;
}

const double ElectrolyteSolution::EffectiveIonicStrength(const VectorXd& n) const
{
	double Ie = 0.5 * (zi.array() * zi.array() * IonicMolalities(n).array()).sum();
	
	return Ie;
}

const double ElectrolyteSolution::StoichiometricIonicStrength(const VectorXd& n) const
{
	double Is = 0.5 * (zi.array() * zi.array() * StoichiometricIonicMolalities(n).array()).sum();
	
	return Is;
}

const VectorXd ElectrolyteSolution::IonicMolalities(const VectorXd& n) const
{
	// The ionic molalities
	VectorXd mi(ions.size());
	
	// The molality factor to trasnform number of moles into molality
	const double molalityFactor = 55.508 / n[iH2O]; 
	
	// Loop over all the ions
	for(unsigned i = 0; i < ions.size(); ++i) mi[i] =  molalityFactor * n[idxions[i]];
	
	return mi;
}

const VectorXd ElectrolyteSolution::StoichiometricIonicMolalities(const VectorXd& n) const
{
	// Check if there is at least one aqueous complex, otherwise return the ionic molalities
	if(complexes.empty()) return IonicMolalities(n);
	
	// The number of ions (Ni) and complexes (Nc) in the electrolyte solution
	const unsigned Ni = ions.size();
	const unsigned Nc = complexes.size();

	// The moles of the ions (ni)
	VectorXd ni(Ni); for(unsigned i = 0; i < Ni; ++i) ni[i] = n[idxions[i]];
	
	// The moles of the complexes (nc)
	VectorXd nc(Nc); for(unsigned c = 0; c < Nc; ++c) nc[c] = n[idxcomplexes[c]];

	// The molality factor to trasnform number of moles into molality
	const double molalityFactor = 55.508 / n[iH2O]; 

	// Return the molality factor times the stoichiometric number of moles of each ion
	return molalityFactor * (ni + vci.transpose() * nc);
}

const VectorXd ElectrolyteSolution::Molalities(const VectorXd& n) const
{
	// The molality factor to trasnform number of moles into molality
	const double molalityFactor = 55.508 / n[iH2O]; 
	
	return molalityFactor * n;
}

const VectorXd ElectricalCharges(const vector<string>& ions)
{
	VectorXd charges(ions.size());
	
	for(unsigned i = 0; i < ions.size(); ++i) 
		charges[i] = ElectricalCharge(ions[i]);
	
	return charges;
}

const MatrixXd StoichiometricMatrixComplexesIons(const vector<string>& complexes, const vector<string>& ions)
{
	// Initialise the stoichiometric matrix of the complexes-ions
	MatrixXd vci = MatrixXd::Zero(complexes.size(), ions.size());
	
	// Get an instance of the database
	const Database& db = Database::GetInstance();
	
	// Loop over all the complexes in the electrolyte solution
	for(unsigned c = 0; c < complexes.size(); ++c)
	{
		// Get the complex info from the database
		const ComplexInfo& complex = db.GetComplexInfo(complexes[c]);
		
		// Loop over all the ions that compose the current complex species
		for(unsigned iter = 0; iter < complex.ions.size(); ++iter)
		{
			// The current ion that compose the current complex species
			const string ion = complex.ions[iter];
			
			// The stoichiometry of the current ion that compose the current complex species
			const double v = complex.stoichiometries[iter];
			
			// Get the local index of the current (ion) in the list of (ions)
			const unsigned i = IndexOf(ion, ions);
			
			// Set the entry (c, i) to the stoichiometry of the current ion that compose the current complex
			vci(c, i) = v;
		}
	}
	
	return vci;
}

//const MatrixXd StoichiometricMatrixComplexesIons(const vector<string>& complexes, const vector<string>& ions)
//{
//	// Initialise the stoichiometric matrix of the complexes-ions
//	MatrixXd vci = MatrixXd::Zero(complexes.size(), ions.size());
//	
//	// The complex map object that will store all the aqueous complexes and their dissociation expression
//	ComplexMap complexMap;
//	
//	// Build the complex map object
//	BuildComplexMap(complexMap);
//	
//	// Loop over all the complexes in the electrolyte solution
//	for(unsigned c = 0; c < complexes.size(); ++c)
//	{
//		// Loop over all pairs (stoichiometry, ion-name)
//		BOOST_FOREACH(auto pair, complexMap[complexes[c]])
//		{
//			// Get the local index of the ion in the current pair (stoichiometry, ion-name)
//			const unsigned i = find(ions.begin(), ions.end(), pair.second) - ions.begin();
//			
//			// Set the entry (c,i) to the stoichiometry of the ion in the current pair (stoichiometry, ion-name)
//			vci(c, i) = pair.first;
//		}
//	}
//	
//	return vci;
//}

